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We study Dirac plasmons and their damping in spatially separated iV-layer graphene structures 
at finite doping and temperatures. The plasmon spectrum consists of one optical excitation with a 
square-root dispersion and TV — 1 acoustical excitations with linear dispersions, which are undamped 
at zero temperature within a triangular energy region outside the electron-hole continuum. For any 
finite number of graphene layers we have found that the energy and weight of the optical plasmon 
increase in the long wavelength limit, respectively, as square-root and linear functions of N. This is in 
agreement with recent experimental findings. With an increase of the number of multilayer acoustical 
plasmon modes, the energy and weight of the upper lying branches also exhibit an enhancement 
with N. This increase is strongest for the uppermost acoustical mode so that its energy can exceed 
at some value of momentum the plasmon energy in an individual graphene sheet. Meanwhile, the 
energy of the low lying acoustical branches decreases weakly with N as compared with the single 
acoustical mode in double-layer graphene structures. Our numerical calculations provide a detailed 
understanding of the overall behavior of the wave vector dependence of the optical and acoustical 
multilayer plasmon modes and show how their dispersion and damping are modified as a function 
of temperature, interlayer spacing, and inlayer carrier density in (un)balanced graphene multilayer 
structures. 

PACS numbers: 72.80.Vp, 63.22.Rc, 71.10.?w, 71.38. ?k 



I. INTRODUCTION 

Graphene^l, a single-atom sheet of graphite, has at- 
tracted considerable attention due to its zero band gap 
electronic structure with massless linear and chiral dis- 
persion of charge carriers^'. These unique characteris- 
tics lead to numerous of extraordinary new phenomena 
in graphene such as the Klein tunneling^, the ambipo- 
lar field effectP^, the finite minimum conductivity^', the 
ultra -high mobility^!, and the anomalo us qua ntum Hall 
effecpESl. With these novel propertie d 11 ' 12 ', graphene 
presents itself as an extremely promisin g ca ndidate- 
material for future generation of electronic d 13 ' 14 '. 

Graphene structures open a new arena also for study- 
ing fundamental many-body interaction phenomena, 
which can play also an important role in graphene based 
electronic devices^. Particularly, electron-phonon in- 
teraction phenomena have become a subject of active 
study in single-layer graphene structures in zercE^HH] 
and finite magnetic fieldd^SHUl Substantial efforts have 
been directed towards the investigation of the linear 
response of d oped g raphene^ and of the charge den- 
sity excitations 25 " 29 and of such complex quasiparticles 
as plasmaron d 30 ' 3 1 ' and plasmon-phonon complexes^!. 
Recently, the experimental realization of graphene 
double-layer structures coupled only via the Coulomb 
interactional 3 ^' has attracted substantial theoretical in- 
terest in studying the double-layer plasmon effects^HUl 
and the frictional dragp^ in two spatially separated 
graphene layerd^M^l ag p 0wer f u i tools for probing inter- 
action effects of massless Dirac fermions. 

One of the key quantities for investigating plasmon 
properties in Fermi liquid theory^ is the dynamical po- 
larizability, which, together with the Coulomb interac- 
tion potential, determines the dielectric function of an 



electronic system and describes its dispersive and dissi- 
pative properties. The recent theoretical calculations of 
the bubble diagram of the Lindhard polarization func- 
tion, Uo(uj, q), in monolayer graphene^^Zl show that 
for an arbitrary bosonic frequency, w, and momentum, 
q, its dependence on the carrier density, n, is weaker, 
IIo(w,g) oc y/n, than the dependence Ho(u}, q) oc n in 
conventional two dimensional electron systems with a 
parabolic energy dispersion. As a direct consequence, 
the frequency of charge density waves, uj p , in a quantum 
system of chiral massless Dirac fermions shows a unique 
carrier density scaling, uj p oc n 1 / 4 , which is in contrast 
with conventional two-dimensional plasma ui p oc n 1 / 2 . 
The weak dependence of the plasmon frequency and its 
weight on the carrier doping level can be a limiting fac- 
tor in the realization of effective control and tunability 
of graphene-based applications, which due to many novel 
properties hold a promise of completely new functionali- 
ties in optoelectronics and terahertz metamateriald^SZI. 

Recently, Ref. 58 demonstrates experimentally a new 
tunable photonic device based on graphene/insulator 
stacks, formed by depositing N — 5 alternating 
wafer-scale graphene sheets and thin insulating lay- 
ers. Measurements^ of light-plasmon coupling in such 
graphene multilayer structures (GMLS) show that dis- 
tributing carriers into several graphene layers enhances 
effectively the plasmon frequency and the magnitude of 
extinction of the light transmission in comparison with 
that in single-layer graphene of the same total density. 
The measured frequency and weight enhancement of the 
in-phase optical plasmon mode follows, respectively, the 
square-root and the linear dependence on the number 
of graphene sheets in GMLS, i.e. \N and 7Y, which 
differs from the weak inlayer density dependence, n 1 ^ 4 
and n 1 / 2 , obtained respectively for the plasmon fre- 
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FIG. 1. (Color online) The plasmon dispersions at finite tem- 
perature T = 0.1T F in GMLS consisting of N = 2,3,4 and 
5 graphene monolayers (shown in bold parallel lines in each 
panel) for the interlayer separation d = 5 nm. The upper solid 
curve represents the optical plasmon mode with the square- 
root dispersion while lower dotted, dashed, dotted-dashed, 
and long-dashed curves represents acoustical branches of the 
plasmon spectrum with linear dispersions. The thin black line 
represents the plasmon mode in an individual graphene layer. 
The doping level in each layer corresponds to the inlayer car- 
rier density m = no = 10 12 cm -2 . The thick dot-dashed lines 
show the boundaries of inter- and intra-subband electron-hole 
continua. All quantities in this figure are dimensionless. 

quency and its weight both in single-I^H^Sl and double- 
layep22JH2J graphene structures. This effect is strikingly 
different from that in conventional two-dimensional sys- 
tems in semiconductors where, if exchange and correla- 
tions effectsP are neglected, the in-phase plasmon mode 
in a symmetrically balanced double-layer system behaves 
similarly to that in an individual layer with the doubled 
density. The observed unusual behavior of the in-phase 
optical plasmon mode in GMLS has been ascribed^ to 
the nonclassical plasmon frequency of massless Dirac 
fermions, which depends explicitly on the Planck con- 
stant h. 

In this paper we present a detailed study of the spec- 
trum of multilayer optical and acoustical plasmon modes 
in GMLS, consisting of a finite number, N, of spatially 
separated graphene layers, where Dirac fermions are cou- 
pled only via many-body Coulomb interaction. Within 
the usual plasmon theory, based on the matrix Dyson 
equation for the Coulomb propagator in GMLS, we cal- 
culate the dispersion relations and Landau damping of 
multilayer plasmons and provide an understanding of 
their dependences on the inlayer electron density and 
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FIG. 2. (Color online) The damping function of multilayer 
plasmon modes at finite temperature T = O.ITf in GMLS 
consisting of N = 2,3,4 and 5 graphene monolayers. The 
curves plotted here corresponds to the plasmon branches rep- 
resented in Fig. [T] with the same formatting. All system pa- 
rameters are taken the same as in Fig. [T] 

the number of graphene layers, N. The plasmon modes 
obtained here for structures with a finite number of 
graphene layers are qualitatively different from the modes 
obtained previously in infinite graphene superlattices^ 
and in structures with a conventional parabolic spectrum 
of electrons®!. We find that the plasmon spectrum con- 
sists of one optical mode with a square-root dispersion 
and N — 1 acoustical modes with linear dispersions in 
the bosonic momentum q and in the long wavelength 
limit all these modes are undamped at zero tempera- 
ture. In this limit we derive analytical formulas for the 
frequencies of multilayer optical and acoustical plasmon 
modes in GMLS. The frequency and the weight of the 
optical plasmons show an enhancement, respectively, as 
y/~N and N, in full agreement with the measurements of 
Ref. [SSI We find that in GMLS the velocity of acoustical 
plasmons at vanishing momenta q is given by a universal 
formula with different effective screening wave vectors, 
which differ from the Thomas-Fermi screening wave vec- 
tor in an individual graphene sheet by numerical factors, 
determined by the total number of graphene layers N and 
by the index of acoustical modes. Furthermore, we have 
presented detailed numerical calculations of the plasmon 
energy spectrum and the damping function in the full 
range of momenta in GMLS, consisting of up to N = 5 
graphene layers. Our comprehensive study includes finite 
temperature effects in different configurations of density 
balanced and unbalanced GMLS. 

The paper is organized as follows. In Sec. II we present 
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FIG. 3. (Color online) The plasmon dispersions in GMLS 
consisting of N = 4 graphene monolayers for four different 
temperatures, T — Tjr. All other parameters are the same as 
in Fig. |I| 



FIG. 4. (Color online) The broadening of the plasmon disper- 
sions. The different curves here correspond to the multilayer 
plasmon branches with the same formatting and the same 
parameters as in Fig. [3] 



the model and use it to describe the multilayer plasmon 
modes in GMLS. Here we derive analytical expressions 
for the plasmon spectrum in TV-layer graphene structures 
in the long wavelength limit. In Sec. Ill, we present our 
numerical calculations for arbitrary momentum and dis- 
cuss the energy dispersions of the optical and acoustical 
plasmon modes and their corresponding damping as a 
function of N. In this section we study the effect of fi- 
nite temperatures on the behavior of multilayer plasmon 
modes in different density balanced and unbalanced con- 
figurations of GMLS. In Sec. IV, we close the paper with 
a summary of our main results. 



II. THEORETICAL MODEL AND 
ANALYTICAL RESULTS 

We consider GMLS consisting of ./V Coulomb coupled 
graphene layers with an equal interlayer spacing d and 
obtain the elementary excitations of the multilayer plas- 
mon modes from the poles of the exact Coulomb Green 
function, V(q,uj). In such GMLS the bare Coulomb in- 
teraction is a symmetric tensor 



= v(q) exp - j\qd) 



(1) 



with respect to the layer indices i, j — 1, 2, . . . , N. Here 
d is the interlayer spacing. The non-diagonal elements 
in Eq. ([I]) represent different interlayer electron-electron 
interactions while the diagonal elements are given by the 
Fourier transform of the intralayer Coulomb potential in 



two dimensional momentum space, v(q) = 2ire 2 /nq. Ac- 
cording to the experimental situation of Ref. [SHI here 
we assume that the dielectric background of the environ- 
ment surrounding graphene layers is homogeneous and 
can be well described by an average effective dielectric 
constant k. The matrix kernel v(q) = Vij(q) determines 
a standard matrix Dyson equation for the exact Coulomb 
Green function 



V(q, u) = v(q) + v(q) ■ fl(q, u) ■ V(q, u) 



(2) 



where Tl(q, u>) is the irreducible polarization function of 
the GMLS. As seen from |2]) the poles of the exact 
Coulomb Green functions are given by the zeros of the 
real part of the scalar screening function 



with 



Si £at(<z, u>) = 



£jv(<?,w) = det ¥ — v(q) ■ fl(q,uj) 



(3) 



(4) 



In experiment the interlayer spacing in GMLS is suf- 
ficiently large so that the electrons in different spa- 
tially separated graphene monolayers are coupled only 
via many-body interlayer Coulomb interaction, i.e. the 
electron tunneling between layers is insignificant and 
we can neglect the non-diagonal elements of the po- 



larizability and take U(q, U)) 



5ijH l (q,uj). In the 



self-consistent random phase approximation^ n^g, to) is 
given in terms of the noninteracting Lindhard polariza- 
tion function U^(q,u>), calculated quantum mechanically 
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from the bubble diagrams in the ith graphene monolayer 
with an electron density rij. This approximation neglects 
vertex corrections^, however, it usually provides a sat- 
isfactory description of the plasmon modes within the 
Fermi liquid formalism. 

In balanced GMLS with an equal electron density in 
all monolayers, = n , the screening function can be 
represented as a product of the form 

N 

e N (q, w) = I] (! - rN(qd)v(q)Ilo(q, uj) (5) 
i=i 

where the set of N different algebraic functions /jy de- 
termine the actual dispersion relations of the plasmon 
modes in GMLS. An important feature of the Coulomb 
interaction is that for an arbitrary number of layers N 
one of these functions /jy in the long wavelength limit 
does not depend on the bosonic momentum q and is given 
namely by the number of layers in GMLS, 

fh(qd) = N . (6) 

Accordingly, determines the frequency of the in-phase 
optical plasmon mode with a square root dispersion with 
q. Moreover, it can be shown that the remaining N — 1 
functions P N with i = 2, . . . , N vanish linearly with q in 
this limit and can be represented as 

f]+\qd) = a l N qd . (7) 

These functions identify the A— 1 out-of-phase acoustical 
plasmon modes with a linear dispersion. The numerical 
coefficients a l N with i = 1, . . . , N — 1 determine the ve- 
locity of the ith acoustical mode. 

As seen from Eq. ^ the effective polarizability for the 
in-phase optical mode is II, n = NU (q, uo). Its imaginary 
part determines the weight of the optical plasmon hence, 
in agreement with the experimental findings of Ref. I58[ 
it is additive in the number A of monolayer graphene 
sheets in GMLS. This conclusion is independent of the 
quantum mechanical properties of the electron motion 
in each graphene layer that determines the properties 
of the Lindhard polarization function n (q, w) of single 
monolayer graphene. In the long wavelength limit the 
in-phase optical plasmon mode can be obtained from the 
high-frequency dynamical (v F q <C oj <C e F ) asymptote 
of the Lindhard function, which can be written in the 
following universal form 

U (q,uj) cx ^-^2 ■ ( 8 ) 

In conventional two-dimensional electron gases with 
parabolic dispersion, the Fermi energy increases linearly 
with the doping level, e F oc hq, therefore the total in- 
phase polarizability given by Eq. ^ is additive in the 
total number of electrons, IL„ cx n to t with n tot — Nuq, 
while the in-phase optical plasmon frequency uj op cx 
■J ntot- The situation is strikingly different in GMLS. 
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FIG. 5. (Color online) The plasmon dispersions at finite tem- 
peratures T = O.ITf in GMLS consisting of N = 2,3,4 and 
5 graphene monolayers (shown in bold parallel lines in each 
panel) for the interlayer separation d — 20 nm. All other 
notations and parameters are the same as in Fig. [l] 

The unique spectrum of Dirac fermions with the con- 
stant velocity implies that the Fermi energy in a graphene 
monolayer, ep = vpkp with kp — ^iruo, is a square 
root function of the inlayer density uq hence the total in- 
phase polarizability of GMLS, again given by Eq. (JsJ) , is 
no longer additive in n to t ■ However in virtue of Eq. (|6| , 
it is still additive in the number of layers N so that we 
have IL„ cx Af^/ij) and u} op cx %/AYio 1 / 4 , which is again 
in agreement with the measurements of Ref. 1581 

Because the dispersion of the acoustical plasmon 
modes is linear, one cannotP^^ use the dynamical 
(fijq — > oo) limit ^ of the Lindhard function in the 
long wavelength limit in order to obtain the plasmon ve- 
locity. Instead, we obtain the dispersion of acoustical 
modes making use of the exact expression^ for the zero 
temperature Lindhard polarization function in graphene. 
Thus, in the long wavelength limit we find the following 
square-root and linear energy dispersions 

i \ I at 9 e2v Fk F i , , l + a l N q TF d 

u p(q) = \ Nq , uj ac (q) = , =«Fg , 

V K yfl + 2a l N q TF d 

(9) 

respectively, for one optical and N — 1 acoustical mul- 
tilayer plasmon modes. Here g = 4 accounts for the 
spin and valley degeneracy in graphene and qxF — 
Ae 2 k F /(n£p) is the Thomas-Fermi screening wave vector 
in graphene. The numerical factors a l N (i = 1, . . . , A— 1) 
of the ith acoustical mode are uniquely determined by the 
total number of graphene layers N in GMLS. For N = 2 



5 




FIG. 6. (Color online) The plasmon dispersions in an unbal- 
anced N = 3-layer graphene structure. The different panels 
correspond to the inlayer carrier density rii — 0, 0.4, 0.8, 1 x 
10 12 cm -2 of the middle layer graphene sheet. All other pa- 
rameters and notations are the same as in Fig. [T] 



FIG. 7. (Color online) The plasmon dispersions in an un- 
balanced N = 5-layer graphene structure with two unroped 
layers. The different panels correspond to the different situa- 
tions of m = ri3 = 0, m = 7U = 0, ri2 = ns = 0, ri2 = rit± — 0. 
All other parameters and notations are the same as in Fig. [T] 



we have only one coefficient a\ = 1 and the above formu- 
las recover the previous results obtained for the plasmon 
dispersions in double-layer graphene structure d 39 * 41 * 42 !. 
For GMLS consisting of N — 3 graphene layers, we find 
for the two coefficients 



1 O 2 

a, = 2 , a 3 



2/3 , 



(10) 



corresponding to the two acoustical plasmon modes. For 
N = 5 the four different velocities of the acoustical plas- 
mon modes in Eq. ^ are given in terms of the following 
coefficients 



1,2 



3,4 



3±V5 



5± 



(11) 



In the limit of small qrFd, the separations between ith 
and jth acoustical modes from the top of the intrasub- 
band electron-hole continuum differ in the long wave- 



Ill. PLASMONS DISPERSIONS IN 
MULTILAYER GRAPHENE STRUCTURES 

In this section we present our numerical calculations 
for the energy dispersions and Landau damping of mul- 
tilayer plasmon modes both in balanced and unbalanced 
GMLS with up to N = 5 graphene layers. We calcu- 
late the energy spectrum from the general formula ^ of 
the dispersion equation with the screening function Q. 
For the diagonal elements of the polarization function 
n o (q,u;) we make use of the exact semianalytical formu- 
las from Ref. for the Lindhard polarization function 
and study also the effect of finite temperature both on 
the plasmon dispersion and damping. The Landau damp- 
ing of the plasmon modes can be well described by the 
broadening function of the respective energy dispersions 



length limit by a factor 
2 

■ ; "* It gives a factor of 73 



("jv/ a 7v) 



(uj l ac - v F q)/ (uj J ac - v F q) = where the function 



9 difference for the 



acoustical branches in GMLS with N = 3 while for = 5 
the difference for the most distant acoustical modes is 
even larger, by a factor of 75 ~ 90. In experiment, how- 
ever, the parameter qTFd is typically not very small and 
this difference is moderate. For the inlayer carrier den- 
sity no = 10 12 cm~ 2 and the interlayer spacing d = 5nm 
taking k = 3.8, we have qTF<i ~ 2 and this gives 72 ~ 3 
while 75 ~ 7. 



3n°(g,w) 
dnu°(q,uj)/duj 



(12) 



(13) 



gives the property of the electron-hole continuum to 
cause damping of the elementary excitations. Thorough 
the calculations we assume that GMLS has a homoge- 
nous dielectric background with the dielectric constant 
k = 3.8. In actual calculations of the dispersion relations 
we set also the imaginary part of n o (q, w) to zercESin the 



G 



screening function Q, which provides a correct asymp- 
tote for the spectral branches of multilayer plasmons in 
the limit of d —> oo. 

In Fig. [I] we plot the dispersion relations of multilayer 
plasmons in graphene structures with iV = 2, 3, 4 and 5 
layers. The plasmon dispersion in a single layer is shown 
also by a thin line for reference. We see a clear shift of 
the in-phase optical plasmon mode towards higher ener- 
gies with an increase of the number of graphene layers 
in GMLS. This enhancement is strongest at wave vectors 
around q = 0.3kp where we find an increase of the dis- 
tance between the multilayer and the single-layer optical 
plasmon modes of about 0.17e p for a double-layer struc- 
ture, N = 2, and of about 0.39sf for GMLS, consisting of 
N = 5 graphene layers. Strong modifications occur also 
for multilayer acoustical plasmons. With an increase of 
N we find an increase of the number of multilayer out-of- 
phase modes and also an enhancement in the energy of 
the upper lying acoustical modes. This enhancement is 
so strong that the top acoustical mode crosses the single 
layer plasmon mode in GMLS for N — 4 and 5 in the tri- 
angular region outside the electron-hole continuum. We 
find that in GMLS consisting of N = 5 (N = 4 and 
N = 3) graphene layers, only 2 (1) from the 4 (3 and 2) 
acoustical branches are suppressed as compared with the 
single acoustical mode in GMLS for N = 2. The behav- 
ior of multilayer plasmon modes in the short wavelength 
limit remains almost unchanged. In this approximation, 
independent of the number of layers, all multilayer plas- 
mon modes for q 3> kp coagulate around the single-layer 
plasmon branch and remain at some distance parallel to 
the boundary of the intra-chirality subband electron-hole 
continuum. 

In Fig. [2] we study the Landau damping of the cor- 
responding multilayer plasmon modes, shown in Fig. [I] 
The energy enhancement, observed for all multilayer 
plasmon modes, accompanies with an enhancement of 
their damping in the number of layers N. It is seen 
that the damping is strong where the plasmon mode 
is far from the boundary of the intra-chirality subband 
electron-hole continuum hence the decay of the in-phase 
optical mode is always larger than that of the out-of- 
phase acoustical modes. This indicates that at finite 
temperatures the finite damping of the multilayer plas- 
mon modes outside the electron-hole continuum is due to 
a leakage of the plasmon weight from the inter-subband 
electron-hole continuum towards low energies and mo- 
menta but from the intra-subband continuum, at which 
boundary, uj — vpq, the damping function vanishes, 
T(q,vpq) = 0. The broadening of a multilayer plas- 
mon mode shows a strong increase at some upturn mo- 
mentum, q c , which is the smallest for the optical mode 
and the largest for the lowest in energy acoustical mode, 
respectively, with the values q c s» Q.2kp and 0.8kp in 
GMLS consisting of N = 5 layers. In the short wave- 
length limit, q 3> kp, the broadening of multilayer plas- 
mon modes, Ti(q), increases approximately linearly in q, 
independent of the plasmon mode index i. 



In Figs. [3] and [4] we study the finite temperature effect 
on the plasmon dispersions and on the broadening of the 
respective plasmon modes in GMLS consisting of N = 4 
layers. It is seen that with an increase of T, the energy 
of all multilayer plasmon modes exhibits a significant 
increase. This enhancement results in a larger broad- 
ening of the corresponding plasmon dispersions because 
the distance of the respective modes from the boundary 
of the intrasubband electron-hole continuum increases. 
The other source of the broadening enhancement is the 
increasing leakage of the plasmon weight with T from 
the intersubband electron-hole continuum. The Landau 
damping increases drastically, especially for the in-phase 
optical and for the top out-of-phase acoustical plasmon 
modes so that the upturn momentum for these modes al- 
most vanishes. It is seen that at temperatures T ~ 0.5Tp, 
the multilayer plasmon dispersions show up to the 15 — 20 
percent broadening in the triangular energy region out- 
side the zero temperature electron-hole continuum. 

The plasmon dispersions in GMLS consisting of up to 
N = 5 graphene layers are shown in Fig. [5] for a larger- 
value of the interlayer spacing, d = 20 nm. One can see a 
clear convergence of the optical and all acoustical multi- 
layer plasmon modes to the single-layer plasmon disper- 
sion. This tendency is so strong that it is hard to distin- 
guish between the multilayer plasmon modes when they 
enter into the inter-subband electron-hole continuum. 

In Figs. 6][7 we study the spectrum of the multi- 
layer plasmon modes in unbalanced GMLS. Fig. [6] shows 
the plasmon dispersions in GMLS consisting of N = 3 
graphene layers for four different values of the inlayer 
density, 712, in the middle layer graphene sheet. As seen 
in the top right panel for the vanishing density ri2 = 
and finite temperatures T = O.lTp, the lowest acous- 
tic mode is slightly separated from the boundary of the 
intra-subband electron-hole continuum. With an increase 
of 7J2 the separation increases along with an overall en- 
hancement of the energy of all multilayer plasmon modes. 
In Fig. [7] we calculate the plasmon dispersions in GMLS 
consisting of N — 5 graphene layers. The dopping level 
in two layers of them is at the Dirac point. It is seen that 
accordingly only two from the four acoustical modes are 
well separated from the boundary of the intra-subband 
electron- hole continuum at low temperatures T = O.lTp. 
Notice that the energy of all multilayer plasmon modes 
depends weakly on the position of the pristine graphene 
layers. In the configuration shown in the bottom left 
panel, the optical plasmon is slightly suppressed while 
the acoustical modes are slightly enhanced in energy as 
compared with the plasmons in the other configurations. 



IV. SUMMARY 

We have theoretically investigated the collective 
plasma oscillations in iV-layer graphene structures where 
Dirac Fermions in different spatially separated layers are 
coupled via Coulomb interaction. We have calculated the 
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energy dispersions and Landau damping of the multilayer 
plasmon excitations by changing the total number of lay- 
ers, the inlayer carrier density, the interlayer spacing, the 
lattice temperature, and the number and positions of un- 
dopped layers in such graphcnc structures. The multi- 
layer plasmon spectrum consists of one optical plasmon 
mode with a square-root dispersion and N — 1 acoustical 
plasmon modes with linear dispersion. We have found 
that in the long wave length limit and for any number of 
graphene layers in GMLS, the energy of the optical plas- 
mon mode and its weight exhibit, respectively, a square- 
root and a linear enhancement with N. This calculated 
behavior is in agreement with the recent experimental 
findings of Ref. 58 : and provide a detailed understanding 
of the observed enhancement effect. The velocity of the 
upper lying acoustical branches of the multilayer plasmon 
spectrum increase with the number of graphene layers. 
The obtained enhancement is largest for the top outer- 



most acoustical mode, which crosses at some value of the 
wave vector, q ~ hp, the plasmon energy of an individual 
graphene sheet. The velocity of the low lying acoustical 
branches exhibits a weak suppression in graphene struc- 
tures consisting of N > 3 layers as compared with the sin- 
gle acoustical mode in double-layer graphene structures. 
The presented numerical calculations provide an under- 
standing of the finite temperature effect on the plasmon 
dispersions, their behavior versus the interlayer separa- 
tion and the inlayer doping level. 
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